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Abstract 



A novel time domain solver of Maxwell's equations in passive (dispersive and ab- 
sorbing) media is proposed. The method is based on the path integral formalism of 
quantum theory and entails the use of (i) the Hamiltonian formalism and (ii) pseu- 
dospectral methods (the fast Fourier transform, in particular) of solving differential 
equations. In contrast to finite differencing schemes, the path integral based algorithm 
has no artificial numerical dispersion (dispersive errors), operates at the Nyquist limit 
(two grid points per shortest wavelength in the wavepacket) and exhibits an expo- 
nential convergence as the grid size increases, which, in turn, should lead to a higher 
accuracy. The Gauss law holds exactly with no extra computational cost. Each time 
step requires 0(iVlog 2 N) elementary operations where N is the grid size. It can also 
be applied to simulations of electromagnetic waves in passive media whose proper- 
ties are time dependent when conventional stationary (scattering matrix) methods are 
inapplicable. The stability and accuracy of the algorithm are investigated in detail. 
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1 Introduction 



In this study a time domain solver of Maxwell's equation in passive (dispersive and absorbing) 
media is developed. The main motivation of this work is to bring methods of computational 
quantum physics into classical electromagnetic theory One of the great advantages of time 
domain methods over stationary (scattering matrix) methods is that a single simulation of 
the scattering of a wide band wave packet can determine basic physical properties of the 
target (e.g., transmission and reflection coefficients) in the entire frequency band covered by 
the initial wave packet. Time domain methods also allow for a unique possibility to observe 
all immediate effects on fields caused by the target or by a surrounding passive medium, 
which greatly facilitates qualitative understanding of the interaction of an electromagnetic 
pulse with media and targets. Another important advantage is that the target geometry 
(or medium physical properties) may vary with time and this time dependence cannot be 
removed by going over to a moving reference frame. Stationary methods are simply inappli- 
cable to these kind of problems. 

From the computational point of view, the proposed approach is based on pseudospectral 
methods. The essential advantages of pseudospectral algorithms over conventional finite ele- 
ment or finite difference schemes in solving differential equations are [1]: (i) the exponential 
versus polynomial rate of convergence as the grid size (or the basis dimension) increases; (ii) 
the absence of dispersive errors and (Hi) efficiency in numerical calculations. Time domain 
algorithms in combination with pseudospectral methods have become the state-of-the-art 
technique in numerical studies of quantum dynamics by solving the corresponding initial 
value problem for the Schrodinger equation (see, e.g., [2]). A typical algorithm entails an 
approximate computation of an object called the path integral (or functional integral) intro- 
duced by Feynman [3]. Here Maxwell's theory in general dispersive media is reformulated 
in the Schrodinger (Hamiltonian) formalism. Then the path integral formalism is applied 
to the initial value problem in Maxwell's theory in passive media to develop a numerical 
algorithm. The main objective of this work is to give a theoretical assessment of the path 
integral based solver of the initial value problem for Maxwell's equations. Numerical tests 
and applications will be discussed elsewhere [4]. 

It is shown here that basic principles of the path integral formalism lead to a true time 
domain algorithm which indeed enjoys the advantages of pseudospectral methods. In par- 
ticular, among the aforementioned features, (i) is provided by the use of the fast Fourier 
method [5] as a part of the algorithm, when applied to media whose parameters do not 
have discontinuities in space, (ii) is a consequence of Nelson's construction [6] of the path 
integral which is embedded in our algorithm, (Hi) is due to the fast Fourier method and 
some analytical results that speed up numerical computations. The algorithm has another 
great advantage over finite difference schemes: The Gauss law is implemented exactly with 
no extra computational cost (Theorem 8.2). For widely used multi-resonant Lorentz models 
of passive media, the algorithm is unitary, meaning that, the energy of a wave packet is 
preserved exactly in dispersive media with no attenuation (Theorem 6.1). It is also uncon- 
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ditionally stable (Theorems 7.1 and 7.3) versus conditionally stable finite element or finite 
difference algorithms [1] (see also [7]). A possible drawback of the algorithm (to be tested 
numerically) is that the use of the fast Fourier method in combination with Nelson's con- 
struction of the path integral might require additional computational costs for boundary 
value problems with complicated boundary geometry. In our approach, conventional bound- 
ary conditions are not imposed on electromagnetic fields. Targets and medium interfaces are 
modeled by discontinuous medium parameters. The problem arises from well known features 
of the Fourier method [5] : Aliasing and low convergence rates for non-smooth functions. In 
this study we offer one possible way to deal with this problem while keeping the Fourier 
basis in the algorithm. Alternative pseudospectral approaches to circumvent the problem 
exist and are mentioned here, but not discussed in detail. 

The basic idea of the path integral approach to solving linear, homogeneous, evolutionary 
differential equations (numerically or analytically) is based on the Hamiltonian formalism. 
In the framework of the Hamiltonian formalism, an original system of differential equations 
is transformed to an equivalent system of first-order (in time) differential equations by ex- 
panding the original configuration space, that is, by going over to a generalized phase space 
where all time derivatives, save for the one of highest order, become independent variables 
[8]. A generic linear homogeneous first-order system can be written in the form 

d t % = n% , **=o = , (i.i) 

where dt stands for the partial derivative with respect to time t, a linear operator 7i is called 
Hamiltonian, while ty t is called a state vector (or wave function). It is an element of the 
generalized phase space of the system and viewed as a collection (column) of the original 
variables and their time derivatives. The generalized phase space is equipped with an inner 
product and becomes a Hilbert space. State vectors are typically vector-valued functions in 
1R 3 , and the Hamiltonian is a differential operator. The choice of the inner product depends 
on the problem at hand. One usually requires that componets of \I/ t are elements of L 2 (M 3 ). 

In general, upon going over to the Hamiltonian formalism, there might occur constraints 
[9, 10] 

C a ^t = , (1.2) 

with C a being a set of linear operators which do not contain time derivatives; a enumerates 
the constraint operators. The constraints must be preserved in the time evolution which 
is described by (1.1). In other words, the solution is sought in the subspace of the Hilbert 
space defined by (1.2). Depending on the type of constraints, there are different ways of 
developing the corresponding path integral formalism. In Maxwell's theory, the constraint is 
the Gauss law, and it is of the "first class" in the Dirac terminology [10]. The characteristic 
feature of a first class constrained system is that 

[H,C a ] ~ C a , [C a ,Cb] ~ C c . (1.3) 

A consequence of (1.3) is that if the initial configuration ^ satisfies the constraints, then so 
does the solution of (1.1). However, after the projection of the Hilbert space spanned by \& t 



3 



onto a finite-dimensional subspace (e.g., a projection on a subspace associated with a finite 
spatial grid as is done in Section 3), which is required for numerical simulations, the involu- 
tion condition (1.3) can be violated causing problems in simulations. For instance, the Gauss 
law is typically violated in any finite differencing approach to simulations of electromagnetic 
wave packet propagation. Special efforts have to be made to ensure the transversality of 
the radiation field in Maxwell's theory, which, in turn, complicates simulation algorithms 
and increases computational costs (e.g., when enforcing the Gauss law in finite difference 
schemes on the grid via the Lagrange multiplier method). It is one of the advantages of the 
proposed path integral based algorithm that the Gauss law can be strictly enforced with no 
additional computational costs for generic passive media (Theorem 8.2). 

The solution to Eq. (1.1) is 

% = exp (tH) * = Ut^o , t>0, (1.4) 

assuming that the exponential of H exists. If the Hamiltonian is time dependent then the 
following replacement has to be made in (1.4) 

exp (tH) -> Texp ^jf drH^j = U t , (1.5) 

where Texp stands for the time-ordered exponential. The operator U t is defined as the 
fundamental solution of (1.1), d t U t = H t U t with U t =o being the identity operator. The fun- 
damental solution has the semigroup property, U tl+t2 = U tl U t2 - The action of the evolution 
operator U t on the initial configuration can be written via its integral kernel, 

%(r) = / dv'U t (v,v')^ (v') . (1.6) 

Using the semigroup property of the evolution operator, the entire time evolution can be 
viewed as consecutive actions of the infinitesimal evolution operator lA& t) where At is a time 
step. If the kernel of the infinitesimal evolution operator is known, then the kernel of the 
evolution operator can be computed as the convolution 

U t (r,r') = dr 1 ---dr n UAt(r,r n )UAt(r n ,r n -i)---UAt(ri,r') (1.7) 

with At(n + 1) = t. The integration variables can be regarded as points = r(tk), where 
tk = kAt, k — 0, 1, ...,n + l, on a path r(r) connecting points r(r — t) — r and r(r = 0) = r'. 
In the limit At — > the convolution (1.7) can be viewed as a sum over all paths connecting 
the initial and final points. This is the gist of the Feynman path integral representation 
of the fundamental solution of (1.1). A nontrivial problem is to find the measure on the 
space of paths. For example, if H = A (the Laplace operator), it can be shown that the 
limit exists, and that the measure coincides with the Wiener measure which has support in 
the space of all continuous, but nowhere differentiable paths (trajectories of the Brownian 
motion) pinned at the end points. In quantum mechanics, the problem is more subtle, but 
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can still be solved [11]. The existence of the proper measure on the space of paths opens 
up an attractive possibility to use Monte-Carlo methods of computing the sum over paths 
which is the gold standard algorithm in quantum and statistical physics. 

However, the present study does not intend to tackle the measure problem for the path 
integral representation of Maxwell's theory but rather offers a solution of a more modest 
problem. Namely, how the conventional way, outlined above, of deriving the path integral 
from the original differential equation can be used to obtain an algorithm for numerical 
simulations of the convolution (1.7) for a small, but finite At. Similar ideas have been applied 
to non-dispersive and/or random media as well as to scattering problems and waveguides [12]. 
Our approach applies to general passive media and goes beyond the eikonal approximation 
of geometric optics and/or the diffraction theory used in earlier works on path integrals 
in electromagnetic theory. The results obtained here are believed to be useful for further 
development of path integral methods in theoretical and numerical studies of propagation of 
electromagnetic wave packets in passive media. 

The idea of numerical simulations follows from (1.7) rather straightforwardly, namely, 

^t+M=U At ^t- (1.8) 

Thus, finding a state of the system in a sequential moment of time amounts to computing the 
action of the exponential of a differential operator Ti on the state at the preceding moment 
of time. Theoretically, it is sufficient to know U^ t up to (At) 2 . The limit At — > in (1.7) 
would not change if we replace the exact infinitesimal evolution operator kernel by such an 
approximation. In numerical simulations, the limit is never achieved. Therefore a higher 
accuracy is required to make errors small. Note that the errors are accumulated as more 
iterations (1.8) are taken. An expansion of exp(At7i) into the power series up to some 
desired order is known to produce unstable algorithms. Yet another obvious drawback is 
the lack of unitarity of the time evolution, that is, if the Hamiltonian is skew-symmetric 
( ant i- Hermit ian, if a complex phase space is used), H* = —H, then U^Ut = 1. In the 
Maxwell theory, as we shall see, the squared norm of ^ t with respect to the L 2 (IR 3 ) scalar 
product (^1,^2) — J dr^l^>2 is proportional to the electromagnetic energy of the system. 
Consequently, for non-absorbing media the unitarity of the time evolution is required in 
simulations to provide the energy conservation. 

We shall apply Nelson's method of obtaining the path integral representation of the fun- 
damental solution of Maxwell's equations for passive media. It is based on the Kato- Trotter 
product formula for the exponential of a sum of two noncommuting operators and the use 
of the Fourier basis to compute exponentials of differential operators. Actually, in practical 
applications, a simpler version, known as the Lie- Trotter product formula, is used (see the 
textbooks [13] for details and references therein). In computational quantum mechanics this 
is also known as the split operator method. It allows one to keep the differential operators in 
the exponential, and thereby, ensures the correct dispersion relation of simulated electromag- 
netic waves. It will be shown that there exists a particular realization of this idea in which 
the Gauss law holds exactly in simulations. In general, the Gauss law can be enforced by the 
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projection operator formalism developed for the path integral representation of constrained 
dynamical systems (for a review see [14] and references therein, a numerical application 
to constrained wave packet propagation can be found in [15]). The idea is to replace the 
Hamiltonian by its projection on the subspace (1.3). If V is the projection operator, that 
is, C a V^ = for any V, V 2 = V and V* = V, then H is replaced by VHV. In Maxwell's 
theory the projection can be implemented in our algorithm with no extra computational 
costs. A significant difference from the quantum mechanical case is that the Hamiltonian 
7i (or its projection) is not normal, that is, it does not commute with its adjoint. This 
feature complicates the stability analysis because the von Neumann criteria is no longer 
sufficient for stability, while still being necessary [16]. Nevertheless, the stability, accuracy 
and convergence analysis of the algorithm can be carried out in rather general settings. 

Since time domain simulations are performed on finite lattices, there is always a moment 
of time when the simulated signal first reaches the lattice boundary. One typically uses 
lattices with periodic boundary conditions. So, the pulse would appear on the other side of 
the lattice interfering with itself, thus leading to totally disastrous results for simulations. 
The problem is usually solved by introducing absorbing boundary conditions (see, e.g., [17] 
(for quantum mechanics) and [18] (for electrodynamics)). It is convenient to set a conducting 
layer at the grid boundary whose conductivity is chosen so that it neither transmits nor 
reflects within the designated accuracy in the frequency domain of the initial pulse. In 
Appendix we briefly describe how such a conducting layer can be obtained. 



2 Maxwell theory in the Hamiltonian formalism 

Dynamics of electromagnetic waves in continuous media is governed by Maxwell's equations 

d t B t = cVxH t (2.1) 
d t B t = -cVxE t , (2.2) 

where c is the speed of light, boldface letters denote three- vector fields in M 3 whose spatial 
arguments are suppressed and the time dependence is indicated by a subscript. No external 
currents and charges (antennas) are included in this study. However, the formalism being 
developed is readily generalized to the case when external time dependent sources are present. 
The electric and magnetic induction vectors, D t and B t , respectively, are subject to the 
constraints (the Gauss law) 

V • D t = V • B t = . (2.3) 

In linear response theory, assumed through out the paper, the electric induction is related 
to the electric field as [19] 

D t = E t + f dr xU Er = E t + P t , (2.4) 

J — oo 
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where xt * s an electric response function of the medium and P t is the medium polarization 
vector. A similar relation can be written for the magnetic field and induction, B t = H t + NL t , 
where magnetization M t is determined by the applied magnetic field and the magnetic 
response function of the medium. 

The relation between inductions and fields must be causal, meaning that the response of 
the medium, P t and M t , can only depend on fields applied to the medium prior to the current 
time t, (e.g., = for t < 0) [19]. A natural way to ensure the causality is to require that 
the response function satisfies a differential equation. In other words, the response function 
is assumed to be the fundamental solution of some time evolution differential equation. 
This differential equation can be obtained from a particular physical model of the medium 
in question. A popular model is the multi-resonant Lorentz model. Let D w and E w be 
the Fourier transforms of the electric induction and field. Then from (2.4) it follows that 
D w = e w E w . The dielectric constant in the Lorentz model has the form 

N uj 2 

^ = l + J2~ ' ( 2 - 5 ) 

and Mf = 0. The physical meaning of the model is transparent. The medium is assumed 
to be made of N sorts of damped harmonic oscillators with resonant frequencies oo a and 
damping coefficients j a . Parameters u pa , called the plasma frequencies, are proportional to 
coupling constants of the oscillators to the external electric field (the electric dipole coupling) 
and also depend on the density of oscillators of the sort a. The density may vary in space. 
So uo pa are assumed to be functions of spatial coordinates. In an empty space, oo pa = 0. If 
the resonant frequency is zero, the one-resonant Lorentz model is equivalent to the Drude 
model of metals [19]. In the Lorentz medium the magnetic response function is zero, while 
the electric response function can easily be found by taking the Fourier transform of (2.5). 
Its explicit form is omitted here because it will not be used. The medium polarization is 
determined by a set of second-order differential equations 

N 

d 2 t V a t + 2 la d t V a t + uft» = u&E, , P t = P* ■ (2-6) 

a=l 

Together with Maxwell's equations, Eq. (2.6) form a system of sought-for causal evolution 
equations which are to be transformed into a system of first order equations by means 
of the Hamiltonian formalism. In finite difference time domain numerical schemes, the 
Hamiltonian formalism for the Lorentz model has been used in [20] to study propagation of 
an electromagnetic pulse in homogeneous Lorentz media. 

In our approach no boundary conditions are imposed on electromagnetic fields at medium 
and/or target interfaces. The latter are modeled by spatially dependent couplings of media to 
electromagnetic fields which are included into the system Hamiltonian. At any interface, the 
couplings are allowed to have discontinuities, or, from a physical point of view, they remain 
smooth but change rapidly, \ w \Vuj p \/uj p >>> 1, at the interface, where \ w is a typical wave 
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length of the incoming wave packet. The conventional boundary conditions are automatically 
generated by the dynamics [19]. Thus, the initial value problem is solved in L 2 (M 3 ) for every 
matter and electromagnetic field component. This implies that the energy of the propagating 
wave packet remains finite (in contrast to the scattering matrix approach based on plane 
wave solutions). 

Let us now formulate the initial value problem for a generic passive medium and then 
apply the formalism to multi-resonant Lorentz models. Combine the fields, inductions and 
medium responses into columns: 

*=(*)■ *■=(&)■ 

Assuming linear response theory, one can write for the Fourier transforms 

4>* = %,i>l , (2.7) 

where the Fourier transform of a general response function, Xun nas to satisfy a dispersion 
relation that ensures causality (like the Kramer-Kronig relations for the dielectric constant) 
[19]. For anisotropic media, \u is a symmetric matrix acting on components of electromag- 
netic fields. With this type of generality all possible media are covered as long as linear 
response theory is valid. The response function Xuj can either be modeled or measured and 
tabulated in some frequency range of interest (determined by the frequency bandwidth of 
the initial wavepacket), say, oo G [c^i,^]- Next, the components of xj 1 are expanded in 
a basis of suitable orthogonal polynomials. An optimal expansion is often achieved in the 
Chebyshev polynomial basis. Chebyshev polynomials are defined in the interval [—1,1] so 
a corresponding rescaling and translation of [a^u^] must be done. By taking the Fourier 
transform of xZ^u = we obtain the desired differential equation 

N 

5>n«f = ^f , (2-8) 
n=0 

where uo p = w p (r) plays the role of the coupling constant between matter and electromagnetic 
fields. The order N is determined by the highest order of polynomials used to approximate 
XZ 1 - The expansion coefficients Xn and the coupling oj p are matrices for anisotropic media. 

The basic idea of the Hamiltonian formalism is to convert the system (2.6) or (2.8) into 
a system of first-order differential equations by introducing auxiliary (matter) fields. The 
number of such fields is determined by the order of the original evolution equation for matter. 
For instance, in the case of the multi-resonant Lorentz model, there are N fields P" each of 
which satisfies a second order differential equation. In the Hamiltonian formalism one would 
have 2N real vector fields, j = 1,2,..., 2N. A simple possibility is to set 

{u pa /u a ) er i , (2.9) 
u a e t a , (2.10) 

-2 la g a - vatf- 1 + u pa E t . (2.11) 



p? = 



2a- 
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The reason of inserting the factor u pa /uj a in the definition (2.9) of the auxiliary fields will 
be evident from what follows. Note that the medium polarization P t must be zero in empty 
space where u pa = 0. The factor in (2.9) simplifies the energy conservation and stability 
analysis. 

For the Lorentz model there is another convenient way to introduce the Hamiltonian for- 
malism by using N complex vector fields Ct which satisfy the first order differential equation 

dtCt = KCt-iUpaKt, (2.12) 

p? = ^ (c+ca , (2.i3) 

where A a = — 7 a + iv a and v a = \Joj\ — 7^ . This representation is defined only if 7 a < uo a 
(i.e., the attenuation is not high). From the numerical point of view, solving a decoupled 
system of N first order differential equation and taking complex conjugation (denoted here 
by an over bar) is less expensive than solving an original system of differential equations to 
compute the medium polarization. 

Returning to the general case, we introduce a set of auxiliary fields £ t to convert (2.8) 
into a first-order system, 

dS = nl^ t + VMF^[ • (2.14) 

The operators and Vmf are determined by the details of going over to the Hamiltonian 
formalism. We shall call 7i M the matter Hamiltonian; it governs time evolution of the 
medium when no external fields are applied. The index F indicates that the electromagnetic 
degrees of freedom are described by fields, not inductions. We shall see shortly that the 
matter Hamiltonian depends on whether ipj: or ip[ is used as independent electromagnetic 
variables. The matrix Vmf describes the coupling of matter to the electromagnetic fields, 
which is emphasized by the index MF (matter-to- field coupling). We introduce a linear time 
independent operator 1Z that acts in the space of auxiliary (matter) fields so that 

^ = n^ t , (2.15) 

that is, the (response) operator 1Z maps a given configuration of auxiliary fields onto the 
corresponding physical response field. It depends on the definition of the matter fields (cf. 
(2.9) and (2.13)). A passive medium is not excited, £ 4 = 0, if no external electromagnetic 
field is applied; that is, the initial condition for Eq. (2.14) is such that it has only the trivial 
solution whenever ipf = 0. Under this condition, the solution of (2.14) reads 

6 = f dre^MVMF^ • 
J — 00 

Hence, the linear response operator in (2.7) is the Fourier transform of the operator 

Xt = 0tKe m MV M F , (2.16) 
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where t is the Heaviside function. Or, vice versa, TZ, 7i M and Vmf must be chosen so that 
the Fourier transform of Xt defined by (2.16) coincides with the known response function Xuj 
of the medium in a designated frequency range. 

Maxwell's equations without external currents can be rewritten in the Hamiltonian form 

defi = n F V t - dti>? = n F V t + v™6 . (2.17) 

The field-to- matter coupling Vfm and the field Hamiltonian Tip are deduced from (2.14) by 
acting on the latter by the operator 7Z, which yields 

V FM = -nn M , (2.18) 

n F = (_j^ x c ^ x ) -nvMF^Ho-nvMF . (2.19) 

It is always possible to set up the Hamiltonian formalism so that TZVmf = and, hence, 
H.f — H.Q- It is not difficult to verify that this holds for the Lorentz model discussed above. 
In the general case, the standard procedure of going over to the Hamiltonian formalism [8], 
where components of £ t are identified with time derivatives of the response field, Q ~ d^ipj 1 , 
leads to the same result that TZVmf = 0. Thus, without loss of generality, the last term in 
the field Hamiltonian (2.19) can be omitted. 

The auxiliary matter and electromagnetic fields (or inductions) are unified into a larger 
column 



The wave function tyf satisfies the Schrodinger equation 

mf = H F *( , H F = . (2.21) 

which has to be solved with the initial field configuration ipf =Q = ipo, while the matter fields 
are assumed to be zero at the initial moment of time, £t=o — 0, e.g., the initial wave packet 
is localized in an empty space region. Equations (2.14) and (2.17) are equivalent to (2.21). 
In a similar fashion, one can derive the Schrodinger equation for Note that 

*t' = S*f , S=(l , S-i= (J , (2.22) 

Hence, 

d^l = H^l , Ti 1 = Sn F S~ 1 . (2.23) 
The corresponding blocks of Ti 1 have the form 

Hi = Ho , V M i = Vmf , (2.24) 

Vim = VFM + nn M -n n = -n n , (2.25) 

ri M = Hm-VmfTZ. (2.26) 



10 



To simplify Vim, Eq. (2.18) has been applied. Observe in (2.26) the aforementioned de- 
pendence of the matter Hamiltonian on the representation of electromagnetic degrees of 
freedom. The use of either (2.21) or (2.23) in numerical simulations has its own advantages 
and disadvantages which are discussed below. 

As an example, we give an explicit form of the Hamiltonian for the Lorentz model when 
the auxiliary field are defined by (2.9) 

Vfm = (Vfmij Vfm2, ■ ■ ■ , Vfmn) , VFMa = ( n n pa J , (2.27) 



x o 

Vmf = ~V* FM , (2.28) 

H F M = diag (H M1 , H M2 , • • • , H MN ) , ft M a= ) ' ( 2 -29) 

where diag indicates that the corresponding matrix is block-diagonal with blocks listed in 
the order from the upper left to lower right corners. Note that the matrices Vfmci and 
H-FMa act on a six-dimensional column (£j a_1 , £ t a )*- Therefore they should be understood 
as composed of 3 x 3 blocks. Each block is obtained by multiplying the unit matrix by the 
number indicated in place of the block in (2.27) and (2.29). 

Our final remark in this section concerns "canonical" transformations in the Hamiltonian 
formalism. As has been pointed out, the auxiliary fields are not uniquely defined. There is 
a freedom of making general complex nonsingular linear transformations such as 

&->Sm&, det<S M ^0. (2.30) 

If the infinitesimal evolution operator U^ 1 ^ = exp(At7i^ jP ' 7 - ) ) is computed with one choice of 
the auxiliary fields, a simple similarity transformation, like the one in (2.23), would allow us 
to compute it in any other basis of auxiliary fields. This is an important observation because 
the auxiliary field basis can be chosen in a way that facilitates computation of the evolution 
operator (e.g., to improve the convergence rate or speed up simulations). For instance, in 
the complex representation (2.13) of the auxiliary fields in the Lorentz model, the matter 
Hamiltonian is diagonal. The corresponding transformation of the auxiliary fields is given 
by 

To transform the whole system into this representation, the Hamiltonian 7i F is replaced by 
S~ 1 H, F S and the wave function tyf by S^ F where S is block-diagonal with the unit matrix 
in the upper left (field) corner and with Sm in the lower right (matter) corner. 



S M (fa . (2.31) 



3 The grid representation of Maxwell's theory 

Consider an equidistantly spaced finite grid with periodic boundary conditions. Let Ar 
be the grid step and n be a vector with integer valued components. Then the dynamical 
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variables are projected onto the grid by taking their values at grid points r = nAr, 



*?(r) -> tt?(nAr) , (3.1) 

where Q denotes the representation, / or F. For simplicity, a cubic grid is assumed here. 
It is straightforward to generalize the discussion to a generic rectangular grid. Consider a 
discrete Fourier transformation associated with the grid [5, 21] 

^ t Q (nA;o) = ^nn^ t Q (n'Ar) , T*T = TT* = 1 , (3.2) 

n' 

where the dual lattice step is ko = 27r/Ar. The grid spatial size L and step must be 
chosen so that the Fourier transform of the initial wavepacket has support within the region 
k € [kmin, kmax] where k = |k|, k max = k and k m i n = 2n/L. The Hamiltonian H Q is split 
into a sum 

H Q = + V Q , (3.3) 

where all the spatial derivatives are included into Hq and V Q contains multiplications by 
position dependent functions. This is always possible for the Hamiltonian described in the 
preceding section. The operator V Q is projected naturally 

V°(r)tf?(r) -> V Q (nAr)^(nAr) . (3.4) 

Consider in the Fourier basis, Hq(V) — > (ik). The projection is then done via the 
discrete Fourier transform 

H Q (V)^(r) _ - V (^) nn , H (m'A;o)^?(n'A;o) . (3.5) 



r=nAr 



In what follows, the rules (3.4) and (3.5) define the action of the operators V Q and and 
their functions on any state vector. The action of a product of V Q and Hq on any state 
vector is understood as consecutive actions of these operators according to the rules (3.4) 
and (3.5), in the order specified in the product. 

The projection (3.5) as well as any action of on state vectors is performed by the 
fast Fourier method [5]. It requires -/Vlog 2 iV elementary operations (flops) with N being the 
grid size. In finite differencing schemes, the action of 7i® on a state vector would require 
mNd operations where the integer m depends on a particular difference scheme used to 
approximate derivatives, and Nj is the grid size used in the differencing scheme. It should 
be noted that, as shown below, the use of the fast Fourier transform eliminates the phase 
error (because the correct electromagnetic dispersion relation is preserved) and operates at 
the Nyquist limit. These two features allows one to reduce substantially the grid size as 
compared with that in a finite differencing scheme, while providing the same accuracy in 
simulations. Recall that, in scattering problems, the phase of the return signal contains the 
most significant information about the target. So, in practice, grids in finite differencing 
schemes are significantly larger (more dense) than grids used in the fast Fourier method 
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in order to reduce the phase errors in the former. Needless to say, the advantage of the 
fast Fourier method in reducing the phase error becomes even more significant in higher 
dimensions because N d /N = (n^/n) D where rid and n are the corresponding numbers of 
grid points per shortest wave length in the wave packet, and D is the grid dimension. The 
Nyquist limit is n = 2, while rid is of order 10 or higher. 



4 The split operator method 

Let denote the L 2 (M 3 ) norm of the wave function, or the Euclidean norm of the corre- 
sponding vector (3.1) in the grid representation. One possible way to compute numerically 
the path integral (1.7) is based on the Kato- Trotter product formula [13] 



lim 



-oc 



for a general \1/ and under certain assumptions about the linear operators A and B in the 
Hilbert space spanned by ^. For our purposes it is sufficient to note that for bounded 
operators, (4.1) always holds and is known as the Lie- Trotter product formula. In the grid 
representation, which would always be assumed, unless stated otherwise, operators A and B 
are finite matrices and, hence, bounded. 

Let us apply (4.1) to the split (3.3), meaning that the operator is used in place of 
A (or B) and, respectively, the operator V Q is used in place of B (or A). The infinitesimal 
evolution operator in (1.8) can be approximated by the first term in the following expansion 

= e M(A + B)=gQ t + At 3 WAt (42) 
gQ t = e AtA/2 e AtB e AtA/2 (43) 

W A , = -±([A,[A,B]]-2[B,[B,A]]) + 0(At) . (4.4) 

By making n larger while keeping nAt = t fixed, the strong convergence in (4.1) guarantees 
that the error can be made arbitrary small for any initial state, 

l(W&) n *o -(O n *o|-0 (4.5) 

as n — > oo for any ^o- The numerical iteration algorithm is then based on the replacement 
of the exact evolution (1.8) by the approximate one 

= ■ (4-6) 

The quantity tAt 2 | W^^ol/I^ol can be used to roughly estimate the accuracy of the al- 
gorithm. A more detailed accuracy analysis is given in Section 8. By making use of the 
Campbell-Hausdorf formula for the exponential of the sum of operators it is possible to ob- 
tain the symmetric product formula in (4.2) to approximate U^t up to any desired order 
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in At, presumably achieving a higher accuracy [22]. This would come at the price of hav- 
ing more exponentials in the symmetric product Q® t . In numerical simulations, one should 
keep in mind that computational costs of decreasing At in the third order split (4.2) (i.e., 
increasing the number of steps in the time evolution) might be less than those of computing 
a lesser number of actions of Q® t in higher order splits. So, the higher order splits are not 
always optimal to achieve a better accuracy [2] . 

On the grid, the action of the amplification operator Q® t is computed according to the 
rules (3.4) and (3.5) applied to, respectively, exp(AtV Q ) and exp(At7i^). Explicit formulas 
for the exponentials of the corresponding operators can be worked out in the field and 
induction representations. If the fields are used as independent variables, then a natural 
choice is 

'Ho \ , / V FM 



The matter Hamiltonian Ti^j can also be transferred into V F if so desired. This rearrange- 
ment affects the accuracy of the method, meaning that the operator (4.4) would change. In 
turn, a rearrangement of operators in the split can be used to improve the accuracy. We 
shall discussed this issue later. Using the Taylor series we infer that 

«p(<Hf) = («»*«•> exp( ° w ,,) , (4.8) 

exp(tfto) = 1 + fcos(ct v /z A) - ll V± + 5111 ( ct v^) Hp ? (4g) 
L J cy— A 

where A = V ■ V is the Laplace operator and V± = 1 — V(A -1 V- ) is the projector 
on transverse fields, that is, V±E = E if V • E = and P±E = if the vector field 
E is conservative, E = V0. The projector V± can be omitted in (4.9) if it is known 
(e.g. from a theoretical analysis of the system) that the fields remain transversal in due 
course. In this case, the two first terms in (4.9) are equal to cos(dV— A). The action of 
exp(t7i ) is computed by the fast Fourier transform according to (3.5). In the Fourier basis, 
—A = k 2 = n 2 ^ 2 ,. Note also that the Fourier transform of the fields ip F is required, while 
the auxiliary fields remain in the grid basis all the time. The exponentials of Ti.^ and V F 
can either be computed analytically for simple models like the Lorentz model, as is shown 
Section 5, or, in general case, by direct diagonalization at each grid site. 

Alternatively, the following approximation of the exponential of an operator can be used 

AtB l + Atg/2 + Agg/12 l + Atfi/2 3) (4 1Q) 

l-At£/2 + At 2 i3 2 /L2 + 1 ' 1 - AtB/2 + [ } ' { ' 

If the matrix B is anti-hermitian, then the approximations (4.10) of the exponential of B 
retain unitarity, which is important for stability of the split algorithm (see Theorems 7.1 
and 7.3). On the other hand, costs of computing the inverse matrices in the right hand side 
of (4.10) can be less than those of computing the exponential. 
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If the inductions are used as independent variables, then a natural choice of the split 
would be 

»'=«.'+ v=(* ~T) + U, «{,)■ ^ 

Making use of the Taylor expansion again we deduce that 

exp(W ') = ( eXp( Q mo) [1 " ex P|^)] U \ . (4.12) 



and, similarly, 



exp(tV 7 ) - ( {n i M) -i [eM l n i M) _ 1]Vmi exp(mL) ) • (413) 



Now we can compare the two splits. The split (4.7) has an advantage over (4.11) because 
it requires less calls of the fast Fourier transform. Indeed, in the former the fast Fourier 
transform is called only for the fields ipf. As it follows from (4.12), the operator exp(tH ) 
acts on both the inductions ty\ and the auxiliary fields. Hence the fast Fourier transform must 
be called for the entire column tyj. . If the number of auxiliary fields is large, there might be 
a substantial difference in the computational speed of two algorithms. The latter, however, 
depends on the choice of the matter fields which, in turn, determines 1Z and therefore the 
number of calls of the fast Fourier transform. Note that the matter fields can always be 
chosen in such a way that only one of the components of £ t specifies the response field 
Thus, the canonical transformation (2.30) can be used to reduce the number of calls of the 
fast Fourier transform. If 1Z is chosen so that it depends on position, multiplication of the 
matter fields by TZ must be done before calling the fast Fourier transform. A significant 
advantage of the split in the induction representation is that the Gauss law can be exactly 
fulfilled without altering the algorithm (see Theorem 8.2). 

In empty space either of the splits reproduces an exact solution of Maxwell's equations 
for any period of time t, provided the initial pulse is bandwidth limited. Indeed, on the 
grid, the initial wave packet is a superposition of a finite number of plane waves. Thanks to 
the linearity of the theory, each Fourier mode is evolved exactly, that is, without any phase 
error, by exp(t7i ) for any t > 0. As final remarks in this section, we note that the algorithm 
can operate at the Nyquist limit: Two grid points per shortest wavelength in the initial 
wave packet [5]. Yet, for the multiresonant Lorentz model, it is unconditionally stable (see 
Theorems 7.1 and 7.3). These features cannot be achieved in any finite difference scheme. 



5 A multi-resonant Lorentz model 

An analytical expression for the exponents of and V Q helps to reduce computational 
costs. Here such analytical expressions are derived for multi-resonant Lorentz models. Let 
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us take first the field representation. Due to the block diagonal structure of TC^ we get 

exp(m^) = diag (exp(tHMi), exp(W^ 2 ), exp(tHM N )) , (5.14) 



exp(tW£j = e" 



-7o* 



_ sinhz/ a t , F x 
cosh v a t H z (n Ma + 7a j 



(5.15) 



where u a = (7^ — u^) 1 / 2 . The exponential (5.15) is easy to compute by expanding H F ta in 
the Pauli matrix basis, which is also a basis for the Lie algebra su(2), and then by using 
the well known formula for the exponential of a linear combination of Pauli matrices. For 
small attenuation, 7 a < oo a , we get i> a = %v a . The hyperbolic functions in (5.15) become 
trigonometric ones and i> a is replaced by v a . The eigenvalues of the matter Hamiltonian are 
A a = — 7 a ±i> a . Hence, Re A a < and amplitudes of the matter fields are always exponentially 
attenuated as t — > 00, unless 7 a = leading to Re A a = 0. 

Computation of exp(tV F ) is a bit more subtle. We start with the observation that the 
characteristic polynomial of V F has a simple form 

N 

det (V F - A) = A 2JV (A 2 + c 2 ) , ul = £ < • (5.16) 

a=l 

This can be proved either by a direct computation or by mathematical induction. So, V F 
has 2N zero eigenvalues and two non-zero ones, A = ±iu p . Let X be the eigenvector of V F 
corresponding to the eigenvalue iu p . Its components have the form 

X i = < K + lu0 M > i = 1,2,..., 2(N + 1) , 

so that X • X — 1 and X • X = X • X = where the dot denotes the Euclidean scalar 
product. The skew-symmetric matrix V F has the following spectral decomposition 

V F = iu p {X ®X - X ®X) . (5.17) 

Taking the square of (5.17) we also infer that 

X <g> X = -co; 2 (V F ) 2 - iu- 1 V F . (5.18) 

The exponential of (5.17) is obtained via the Taylor series and making use of (5.18). The 
final result reads 

exp(tV) = 1 + ^ V F + 2 f ™Mg) V (v *y . (5 . 19) 



In the induction representation, an explicit formula for exp(t7i^) is not that simple. 
To avoid unnecessary technicalities, we limit the discussion to the simplest case of the one- 
resonant Lorentz model. We choose the matter fields so that $,\ = P t and £ 2 = d t £\. In this 
case, non-zero elements of the matter Hamiltonian are 7Ym 12 = 1, T^ira = ~ UJ o ~~ u pi an< ^ 
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7~C I M2 2 = ~~ ^7 • The coupling matrix Vmi has only one non-zero element, Vmiix — w 2 Here 
ujq is the resonant frequency, 7 is the attenuation constant and uj p is the plasma frequency. 
Using the Pauli matrix basis again, we find that the expression (5.15) holds for exp(tH T M ) if 

we replace in it v a by v v = ^joo 2 + oo 2 — j 2 , 7 a by 7 and Ti^ by Tt T M . The lower left corner 

of (4.13) has the form 

(^/^(exp^)-!^/" (I °) • (5.20) 

A further simplification can be achieved by going over to the complex representation (2.12) 
of the auxiliary fields in which the matter Hamiltonian is diagonal. The transformation rule 
is explained in the paragraph after Eq. (2.31). 



6 Energy and norm conservation 

Consider the JL 2 (1R 3 ) norm of *?, \^f\ 2 = J dr^f*^f = (#^,* t Q ). In the grid representa- 
tion, the norm coincides with the corresponding (complex) Euclidean norm, fiff , tyf ) = Ar 
J2 n l I / ^*(nAr)\l/^(nAr) where the sum is taken over all grid sites. By taking the time deriva- 
tive of \^ff\ 2 and using the evolution equation (1.1), it is not hard to deduce that the norm 
is conserved, provided the Hamiltonian is anti-Hermitian 

For a generic passive media this is not the case. So the norm is generally not conserved 
in contrast to the quantum mechanical case. However, we shall see that in the case when 
the matter evolution is described by second order differential equations in time and no 
attenuation is present, the norm coincides with the system energy and is conserved. In 
numerical simulations, this important property can be used to help to control the accuracy. 

Consider multi-resonant Lorentz models with no attenuation, 7 a = 0. We start with 
the observation that the field and matter evolution equations can be obtained from the 
variational principle for the action 



JdtL = Jdtjdv \ (E 2 - B 2 ) + \ £ ({d t V a t f - u%0?) + £ u pa V a t ■ E, 



, (6.2) 



where the electromagnetic degrees of freedom are described by vector and scalar potentials, 
respectively, A t and tp t , so that E t = — V</?j — d t A t and B ( = Vx A t . The units are chosen 
in this Section so that c — 1. The polarization of the medium is expressed via the matter 
fields as P t = J^a^pa^t- The least action principle for the scalar potential (p t leads to the 
Gauss law, V-D t = 0, for the vector potential A t to the Maxwell's equation, d t T>t = V x B t , 
and for the matter field to the medium polarization evolution equation of the Lorentz 
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model with no attenuation, 7 a = 0. The second Maxwell's equation and the Gauss law for 
the magnetic field follow from the relation B ( = V x A ( by taking its time derivative and 
divergence, respectively. The energy of the system coincides with the canonical Hamiltonian 
which is obtained by a Legendre transformation [8] of the Lagrangian L for the velocities d t A t 
and <9 t #". Doing the Legendre transformation, we find the canonical Hamiltonian (energy) 
of the system 



Et = \ I <lr 



E t 2 + B? + £ « 2 + c^tff ) 



\\*?\\ (6-3) 



where 7r" = 5L/5(d t 'd a ) are canonical momenta of the matter fields. To get the last equality 
in (6.3), we have used the relations £ 2a = ir® and $, 2a ^ 1 = cj a $" which follow from comparison 
of the canonical Hamiltonian equations of motion for the canonically conjugate variables 
and 7r" and Eqs. (2.10) and (2.11) with 7 a = 0. Note that the canonical momentum 
conjugate to the vector potential A t coincides with — D t = — E t — P t , not — E t in this 
system. Therefore, the coupling between the electromagnetic and matter degrees of freedom 
is included into the term E 2 = (D t — Pt) 2 of the canonical Hamiltonian. Equation (6.3) 
becomes the conventional expression for the electromagnetic energy in a passive medium 
[19] when 7r" and are replaced by the corresponding solutions of the equations of motion 
with initial conditions 7Tq = i9q = 0. The energy conservation can be deduced either from 
the Noether theorem (because E t is the Noether integral of motion corresponding to the 
time translational symmetry of the action) or directly from the norm conservation of V&f 
(because the evolution operator exp(tH F ) is unitary when j a = 0). 

In numerical simulations, an exact unitary evolution operator U2 t is replaced by its 
approximation Q® t . However, the energy remains conservative: 

Theorem 6.1. The split algorithm is unitary for multiresonant Lorentz models with no 
attenuation, that is, the split algorithm preserves the energy E t+ At = E t . 

Proof. In the field representation, 7i F * = —Ti. F and V F * = —V F and, therefore, G^t ^ s 
unitary. As a result, the algorithm preserves the initial wave packet energy and the norm, 

I SKI = Km\ = |*f| • (6-4) 
In the induction representation, the energy coincides with the norm of tyj. in the measure 
space. The measure is determined by the transformation law = S^f, 

Et = \^ F ^ F ) = \{^^{) = \W\l ' ^S^S- 1 . (6.5) 

Since H 1 is similar to H F , the Hamiltonian Ti 1 is ant i- Hermit ian relative to the \x scalar 
product, 

n 1 *^ = -fiH 1 . (6.6) 

The norm conservation (unitarity) in the split algorithm requires in addition that the am- 
plification matrix Q^ t satisfies the unitarity condition 

GmVGL = V ■ (6-7) 
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This is indeed the case. To prove (6.7), we show that Ti^ and V 1 satisfy the condition (6.6) 
and, hence, the product of their exponentials is a unitary operator relative to the \i scalar 
product. Consider H 1 = S^S' 1 =H I Q + V I so that H F = S^H'qS + 5" 1 V / 5. For the 
Lorentz model, 

5-^^; °^ = - (s-'nisy . (e.s) 

Therefore Hq satisfies (6.6). From the anti-Hermiticity of H F and (6.8) it follows that 

(s-Ws)* = -s-Ws . (6.9) 

Hence, V 1 also satisfies (6.6). Thus, 

which completes the proof. 

The norm (energy) conservation can be used to control numerical convergence, especially 
when the aliasing problem in the fast Fourier transform is present, i.e., when parameters 
of the medium are discontinuous functions in space. In a properly designed algorithm the 
loss of energy (norm) due to attenuation should be controlled by the symmetric part of the 
Hamiltonian operator 

d t E t = -5>«l*?f = \ V?*?) < , (6.11) 

a 

where V®* = V® = (H Q * + H Q )/2 < (a negative semidefinite operator) which is, in this 
diagonal matrix with nonpositive elements. 



7 Stability of the algorithm 

The norm of an operator 7i is defined as 

\\u\\ = sup \mi\ . (7.i) 
1*1=1 

If the operator is normal, that is, it commutes with its adjoint, then its norm coincides 
with its spectral radius p(Ti), the supremum absolute value of its eigenvalues. In general, 
p(7~0 < ll^ll- A family of amplification operators (matrices) GAt( a ) is called conditionally 
stable if there exists a constant C(r, T) such that [16] 

\WUa)\\<C(r,T), (7.2) 

for all At G (0,r), all < nAt < T for some positive r and T, and all parameters a. 
The unconditional stability of £7 At (a) means that (7.2) holds uniformly in n > for any 
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At > and for all a, that is, C is independent of T and r. Parameters a can be all wave 
vectors k supported by the grid or simply grid values of the position vector x. They can 
also include parameters of the medium. Note that if G& t is not normal, then p(GAt) < II^aJ 
and, therefore, the von Neumann condition p{Gte) < 1 is no longer sufficient for stability, 
while still being necessary. 

Theorem 7.1. For multiresonant Lorentz models, the split algorithm in the field repre- 
sentation is unconditionally stable. 

Proof. We shall prove that 

WQLW < 1 , (7-3) 

which leads to the theorem statement 

\\{qly\\<\\qlt<i- (7.4) 

By definition and making use of the inequality, \\AB\\ < \\A\\ ||£>||, we get 

\\Glt\\ = \\e Mn «/ 2 e AtvF e Am ^ 2 \\ 

< ||e A *</ 2 || 2 (7.5) 

because e AtvF is a unitary operator, so its norm equals 1. The operator e m ° is block-diagonal 
(see (4.8) and (5.14)). The norm of a block-diagonal operator is the maximal norm of its 
blocks. The upper left corner block is given by the unitary operator e m ° whose norm equals 
1. We have then 

\\e m °\\ = maxjl, || e *^«||} . (7.6) 

The norm of the exponential of 7if fa can be found by direct calculation using the fact 
that ||A.|| 2 = 1 1 -AM. 1 1 = p(A*A) and the explicit form of e m M* given in (5.15). For small 
attenuation, oj 2 — 7 2 = u 2 > 0, we define z a = (j a / u a ) sm(u a t) so that the largest eigenvalue 
has the form 

\\ e tn Ma \\2 = \\( e tH Ma y e tH Ma U 



(l + 2z 2 a + 2z a ^T4) ■ (7-7) 



Since z a < 7„t = y for t > 0, the function (7.7) is bounded from above by f(y) = e~ 2y (l + 
2y 2 + 2y^J\ + y 2 ). It is not hard to verify that the derivative f'(y) is negative for all y > 0, 
and that /(0) = 1. Hence, replacing t by At/2, we conclude that 

\\e Am MJ 2 \\ <1 , (7.8) 

from which (7.3) follows. For large attenuation (like in Drude metals), u 2 — 7 2 = — 1/° < 0, 
in (7.7) we get z a = i^a/^a) sinh(i/ a i) = z a (t). For t > the latter relation defines the inverse 
function t = t(z a ). Once again, the derivative of (7.7) with respect to z a can be shown to be 
negative for all positive z a while at z a = the function equals 1. So inequalities (7.8) and 
(7.3) hold in this case too. This completes the proof. 



20 



The proof of Theorem 7.1 given above is not the most economical. However, the idea 
of estimating the norm of the exponential of the matter Hamiltonian in order to investigate 
stability of the algorithm can be applied numerically to systems more general than the 
Lorentz model because is local on the grid, that is, it does not contain derivatives. So 
the exponentials of Tt^ and its adjoint are not expensive to calculate numerically for some 
trial values of At to see if (7.8) holds. 

We give an alternative proof of the unconditional stability in the case of the induction 
representation of the multi-resonant Lorentz model where an analytical expression of the 
exponent of the matter Hamiltonian is too hard to find, not to mention its norm. We shall 
make use of the following obvious lemma. 

Lemma 7.2. Let a vector ip t , t > 0, be a solution of the linear equation dtipt = (H + V)ipt 
where the linear operators 7i and V satisfy the conditions 7i* = ~7i and V* = V < 
(negative semi definite) . Then \ip t \ < \ipo\ for alH > 0. 

The proof follows from an obvious relation 

Wt| 2 = 2(^,V^)<0. (7.9) 

As a consequence we also get 

||e* (H+v) || <1 , (7.10) 

for all t > 0. 

Theorem 7.3. For multiresonant Lorentz models, the split algorithm in the induction 
representation is unconditionally stable. 

Proof. If the attenuation is absent, the amplification matrix Q^ t is unitary with respect 
to the energy scalar product (yU-scalar product) as is shown in (6.10). Hence, ||({?At)"1U = 1- 
When the attenuation is switched on, the unitarity of the amplification matrix can get 
violated only through exp(AtV / ) because the operator exp^AtO.^) remains unitary with 
respect to the //-scalar product. The idea is to prove the unconditional stability with respect 
to the //-norm. The theorem statement would follow from the equivalence of the Euclidean 
and \i- norms. Recall that two norms and |^|' are equivalent if there exist two positive 
constants Ci )2 such that 

Ci|*| < I^T < C 2 \% , (7.11) 

for all ^. All topological properties of the space spanned by ^ are the same in one norm as 
in the other; in particular convergence of a sequence, boundedness of a set, boundedness of a 
linear operator, and uniform boundedness of a family of operators are all invariant concepts 
under a change of one norm to the other. If \^\' = then ||^|| M = ||<S _1 AS|| and 

= is -1 *! < us -1 !! 1^1 , 

= ISttL < ||S|U*L= ||5|| IttL , 
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so that the two norms are indeed equivalent 



H^ir 1 !^! < < \v\ . (7.12) 

Since ||^4W|| = ||^4|| for any bounded operator A and a unitary operator U, we infer that 

iK&nu < wslk = \\ eAtvI \\» = iie At5 " lv/5 ir • (7.13) 

When 7 a = (no attenuation), the operator S~ 1 V I S is skew-symmetric (cf. (6.9)). When 
7 a 7^ 0, the operator iS _1 V / iS acquires an addition which is a diagonal operator with nonpos- 
itive elements as follows from (2.26) and (2.29). Therefore the inequality (7.10) must hold 
for it as a consequence of Lemma 7.2, 

||e A * 5 " lv/5 || < 1 , 

from which the uniform boundedness of the family (Gj± t ) n with respect to the /z-norm im- 
mediately follows. By the equivalence of the two norms (7.12), the family (£?£ t ) n is also 
uniformly bounded in the Euclidean norm, 

MGitU < \\s\\ ws- 1 ]] , (7.i4) 

which completes the proof. 

Comment. The same idea of making use of the norm equivalence, which actually goes 
in line with the Kreiss matrix theorem (its last part) [23, 16], can be applied to analyze the 
stability of the split algorithm for generic passive media. It is not hard to find a quadratic 
Lagrangian local in time such that the corresponding Euler equations describe propagation 
of an electromagnetic pulse in generic non-absorbing media. Due to time translation symme- 
try, the system should have a conserved quantity according to the Noether theorem [8] . This 
integral of motion coincides with the canonical Hamiltonian which is a quadratic form of tyf 
if the linear response approximation is valid. By analogy with the /i-norm, one could try 
to identify the canonical Hamiltonian with the new norm of tyf which is conserved by con- 
struction and, hence, in an attenuation-free medium the corresponding evolution operator 
is unitary. Thus, it would always be possible to arrange the split so that the amplifica- 
tion operator is unitary too. From the physical point of view, it is then naturally expected 
that, when absorption is added to the system, the attenuation operator V® would generally 
satisfy the condition (6.11) because Fourier amplitudes of fields are exponentially attenu- 
ated in passive media. The latter would make it possible to apply Lemma 7.2 to prove the 
unconditional stability of the amplification operator with respect to the norm defined by 
the canonical Hamiltonian along the lines similar to the proof of Theorem 7.3. An obsta- 
cle for this rather natural idea to generalize Theorem 7.3 to generic passive media is that 
the canonical Hamiltonian is not, in general, positive definite. It becomes positive only on 
solutions of the equations of motion for matter fields, which is a rather common feature of 
Lagrangian systems with higher order time derivatives. Thus, the canonical Hamiltonian 
does not always define a positive definite quadratic form in the Hilbert space for a generic 



22 



passive media, and, hence, cannot serve as a new (conservative) norm. The study of con- 
ditions on attenuation-free media under which a positive definite and conserved quadratic 
form does exist goes beyond the scope of this paper since it would require the canonical 
formalism and the Noether theorem for theories with higher-order time derivatives, which 
is rather involved for generic passive media. The question can be addressed more easily for 
each particular medium model of interest. However, the unconditional stability might be 
excessive as far as practical needs are concerned. It is more important to make the split 
algorithm convergent for a generic passive medium. Then one should use the equivalence of 
(conditional) stability and convergence according to the fundamental convergence theorem 
due to Kantorovich [24, 16]. 

Our findings in this latter approach are summarized in the following theorem. 

Theorem 7.4. Suppose that the medium response function satisfies the causality con- 
ditions (that is, its Fourier transform has poles only in the lower half of the frequency plane, 
Imw < 0). Let U t be an exact evolution matrix in the grid representation (as defined in 
Section 3), and Q^t be an amplification matrix in some third order split algorithm. Then 
for band-width limited wave packets the split algorithm is (conditionally) stable and for all 
< n < N, T = NAt, and < At < r, there exist a constant C m , which depends only on 
the medium parameters, and a constant W m , which depends also on r, such that 

WGltW < C m + 5(At,T), (7.15) 
5(At,T) = C m (e Wm ™ 2 - l) = 0(At 2 ) , (7.16) 
\\lUt-Gl t \\ < S(At,T) . (7.17) 

Remark. Before proving the theorem, let us discuss its significance for practical ap- 
plications. Inequality (7.15) implies conditional stability, while (7.17) establishes a relation 
between the accuracy (and convergence) of the split approximation and the uniform bound 
in the stability condition (7.15). By making the time step At smaller, any desired accuracy 
can be achieved during the total (fixed) simulation time T. The latter implies, of course, 
that the grid is assumed to be chosen fine enough (in accord with the Shannon sampling 
theorem) to accurately reproduce the initial pulse configuration via the fast Fourier method. 
Indeed, let = Glt^a be a simulated solution, and ^ t = Ut^Q be an exact solution, then 
from (7.17) it follows that 

|*t - #H < S(At,T)\^ \ = 0{At 2 ) , (7.18) 

for all < t < T and any fixed total simulation time T which is roughly 2L/c where L is the 
simulation box size and c the speed of light. Now we turn to the proof. 

Proof. In the grid Fourier basis, U t ^Q = ^ k ^(k)e* k ' x , where k spans the dual lattice. 
By construction of the Hamiltonian, each Fourier mode ^(k) evolves exactly as in the con- 
tinuum case. Since the medium response function satisfies the causality conditions, Fourier 
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amplitudes of the electromagnetic and response fields as well as their time derivatives are 
bounded functions of time. The amplitudes cannot grow infinitely large because of dissipa- 
tion [19]. The number of Fourier modes is finite on the grid (only bandwidth limited initial 
wave packets are considered) and, hence, |\^| < Cq for all t > because components of the 
auxiliary field £j are linear combinations of the response field and its time derivatives. The 
latter inequality is equivalent to the evolution matrix being uniformly bounded for all t > 0, 

\\U t \\<C m . (7.19) 

Let U At — G&t = At 3 WAi and Wa< = W + O(At) for small At according to a third order 
split (cf. (4.2) - (4.4)). Let W rn = C m sup At ||Wa*|| for < At < t and some positive finite 
t. Using the semigroup property U\ t = UkAt and (7.19) we infer that 



WQmW = Wit - (uit - $lt)\\ 

< WltW + WUlt-GltW (7-20) 

< C m + \\UZ t - (U At - At 3 W At ) n \\ (7.21) 

= C m + \\- At 3 ( X> At( „- fe -i)>VAtZ4A t j + • • • || (7.22) 

\fe=0 / 

< C m + C m [(l + At 3 W m ) n -l)] (7.23) 

< C m + 5(At,T) . (7.24) 



Inequality (7.17) readily follows from comparing the right hand side of (7.20) with those of 
(7.21)-(7.24). This completes the proof. 



8 Convergence and accuracy analysis 



To estimate the accuracy of the algorithm at a fixed finite grid size N, consider the following 
quantity 

(8.1) 



0n(N, At) = I ((««)» - {G%Y) <| / |* | < \{U Q M ) n - {Q. 



<l) n 



which specifies a deviation of the approximate solution from the exact one relative to a 
given norm. Here U^ t is an exact evolution operator. The accuracy estimate fl n (N,At) is a 
norm dependent quantity. The choice of norm is usually determined by practical needs. We 
use the norm related to the electromagnetic energy of the system and investigate, first, the 
behavior of (3 n {N, At) as At goes to zero, while Atn = t remains fixed and does not exceed 
some positive constant, t <T. 

Theorem 8.1. For multi-resonant Lorentz models, there exists a positive constant W Q 
such that 



< At 2 TW Q c 2 n 



1.2) 



where cf — 1 and cj = ||«S|| ||«S 1 ||, for all At e (0, r) and nAt < T. 
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Proof. In the field representation Q = F, ||(W At ) n || < 1 and ||^ t ) n || < 1 for any 
integer n, as a consequence of Lemma 7.2 for the multi-resonant Lorentz model. The same 
inequalities hold in the induction representation if the norm is replaced by the /z-norm. 
According to the split algorithm (4.2)-(4.4), - g% = At 3 W&. Let W Q = sup At \\w2 t \\ 
for At e (0, r) for some positive r (a maximal time step used in simulations). We then have 
the following chain of inequalities that lead to the theorem statement 



ul t ) n -{Gl) n \ 



(8.3) 



n-l 

fc=l 

< (n-l)At 3 ||Wfj (8.4) 

< At 2 TH/ F , (8.5) 

In the case of the induction representation, inequality (8.4) holds relative to the /z-norm. The 
theorem statement (8.2) follows from the norm equivalence (7.12), cJ^I-AH < \\A\\^ < ci\\A\\ 
for any operator A. The proof is complete. 

Remark. In simulations, the continuum limit iV — > oo is never achieved. Hence the 
operators in the split algorithm (4.1) remain bounded versus the unbounded case of (4.1). 
It is known that the convergence rate of f3 n (oc, At) as At — > estimated by the operator 
norm as in the right hand side of (8.1) is no longer of order 0(At 2 ) but rather of 0(V~At) 
(see, e.g., [25] and references therein). For unbounded operators, the estimate (8.5) is not 
valid. This suggests that the convergence rate l3 n (N, At) may depend, even significantly, on 
the initial vector \l/o as N increases. 

In a general case, the quantity 5 (At, T) in Theorem 7.4 determines the accuracy of the 
split algorithm with respect to the norm (7.1) on a finite grid. To make simulation errors 
small, it is sufficient to require that 

\\U At -gZ\\=At 3 \\wZ\\«l. (8.6) 

Making use of (4.4) and the fact that the norm of a matrix does not exceed the maximal 
norm of its blocks, we infer for a multi-resonant Lorentz model that, in order for (8.6) to 
hold, the following inequalities are sufficient: 

U pa At « 1 , UJmaxAt « 1 , UJ a At « 1 , 7a At « 1 , (8.7) 

and, yet another one, 



cAt « 1 



Here oo rnax is the maximal frequency of the initial wave packet. The right hand side of (8.6) 
is a sum of two types of terms. There are terms of the cubic order in numbers (8.7) as well 
as a term linear in (8.8) with the coefficient being quadratic in (8.7). The ratio in (8.8) can 
roughly be estimated from |Vcj pa | < u) pa / Ar with Ar being the grid step. The condition 
(8.8) implies then that the distance traveled by the wave packet during one time step should 
be much smaller than the grid step. 
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To complete the discussion, one should also analyze the accuracy of the Gauss law (2.3). 
Note that the constraints are automatically fulfilled in the continuum theory due to the 
Dirac involution relations (1.3). By projecting the continuum theory onto a finite grid 
and replacing the exact evolution operator by its approximation in the split algorithm, the 
involution relations might be violated, thus leading to errors and potential instabilities of 
the algorithm. A good example of this kind is numerical general relativity (although the 
nonlinearity of the equations of motion plays the central role in generating instabilities due 
to the violation of the Dirac involution relations). 

It is not hard to be convinced that the Gauss law (2.3) is equivalent to the following 
constraint on state vectors 

0^1 = 0, C'=(£ °) , C=(J ^P„, V {l = l-V ± . (8.9) 

The operator V\\ projects a vector field onto its longitudinal component. In other words, 
it acts as the identity operator if the vector field is conservative, and it annihilates any 
rotational vector field (which is the curl of another vector field) . In the field representation 
we get C F = «S _1 C 7 «S with S defined in (2.22). On the grid, the action of the operator 
C Q is defined by the rule (3.5), that is, by (8.9) we understand {TC Q T*)T^ = 0. Thus, 
the Gauss law requires that the Fourier transform of the inductions should not acquire 
components parallel to wave vectors of the dual grid. This is obviously guaranteed if the 
exact evolution operator, 1A® = exp(tH Q ), is used to generate the time evolution because 

C Q H Q = 0, (8.10) 

and, hence, C Q U?^ = C Q ^ = 0. A problem may arise when the approximate evolution 
operator, (G® t ) n , is used to evolve the initial wave packet \E^. From linearity of the system, 
it is natural to expect that the Gauss law violation should be of the same order as the 
accuracy of a numerical solution of dynamical Maxwell's equations. However, we shall take 
a closer look at the problem and find a pleasant result important in practice, which is stated 
in the following theorem. 

Theorem 8.2. Assuming linear response theory for any passive medium, the Gauss law 
holds exactly in the split algorithm in the induction representation. 

Proof. In the induction representation, identity (8.10) is equivalent to two identities for 
the blocks of C®H Q , namely, CHq = and CVjm = 0. The first one is obvious. The second 
one follows from (2.25) established for any passive medium. The key observation is that the 
identity 

C 7 ^ = (8.11) 

holds thanks to the two above identities and (4.11). Indeed, in the Fourier basis (8.11) is 
equivalent to the vanishing of the triple vector product k • (k x A) for some A regular at 
k = 0. Then from (8.9) and (8.11) it follows that 

CV = C\U l -Hi) = . (8.12) 
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As a consequence of (8.11) and (8.12), we infer that 

C 1 {QitfK = ^e Atw <V 2 e Atv/ e Aw o/2 (gi^ & Q = C l [Q 1 ^' 1 ¥ = C'tf = , (8.13) 
which is the statement of the theorem. 

In the field representation the Gauss law can be enforced by means of the projection 
formalism discussed in Section 1. The projection operator is, obviously, V — 1 — C F . Its 
action is computed in the grid representation by the fast Fourier method according to the 
rule (3.5). Without the use of the projection formalism, the accuracy of the Gauss law is 
stated in the following technical proposition. 

Proposition 8.3. Let W = ||W&|| and W c = || [C F , W&] ||C m (l + 6 (At, T)) where C m 
and S(At,T) are defined in Theorem 7.4, then 

\C F {g F At ) n <| /«| < TW c At 2 + At A WW c e TWM \At 2 + T 2 /2) = 0(At 2 ) , (8.14) 

for all < n < N, T = NAt and any positive At. 

Proof. Since C F U F = C F , assuming that the initial state ^f F satisfies the Gauss law we 

get 

c F {Gl t ) n ^ = -c F {(u F t T-(g F t T}v F 

n—l 

= -c F J2(uL) k Kt-G F A t )(&- k < 



k=l 

n—l n—l 



= -At*[C F ,WZ t ]J2(GL) n - k * F + At'Wl t Y,C F {Glt) n - k < -(8-15) 

k=l k=l 

Denoting the left hand side of (8.14) by a n , we infer from (8.15), by taking the norm of both 
sides, that 

n—l n—l 

a n < At 3 \\[C F ,W At ]\\ \\(&- k \\ + At 3 \\W F t \\Y,* n -k , 

k=l k=l 

for n > 1 and ct\ < WcAt 3 . By Theorem 7.4, powers of the amplification matrix Qj^ t are 
bounded. Hence the following recursion inequality holds 

a n < in - l)At 3 W c + At 3 W(a n ^ + a n _ 2 + • • • + a x ) . (8.16) 

Iterating (8.16) n — l times, we deduce that 

a n < (n - l)At 3 W c + At 3 W {(n - 2) At 3 W c + (1 + At 3 W)(a n _ 2 + a n _ 3 + ---a 1 )} 

{n-2 
J2( n - 2 - fc)(l + At 3 W) k + (1 + At 3 W) n ~ 2 
k=0 
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One can find an explicit form for the sum in the latter equation. However, it is a cumbersome 
expression. For practical purposes, we give a simpler estimate which is stated in (8.14). First, 
factor out (l + At 3 W) n ~ 2 in the brackets, and then use obvious inequalities (l + At 3 W)~ k < 1 
and (1 + At 3 W) n ~ 2 < exp(TTVAt 2 ), which leads to (8.14). 

In the case of the Lorentz model, Wc = \\ [C F , W^J || because all powers of the amplifica- 
tion matrix are uniformly bounded by 1. For small At, a good estimate can be obtained by 
computing W c for At = using (4.4). 

The convergence rate as the number of grid points iV increases is determined by the con- 
vergence rate of the fast Fourier transform which is exponential versus polynomial in finite 
difference schemes, provided parameters of the medium are smooth functions of position 
[1, 5]. As is well known from Fourier analysis, the convergence rate can be affected for func- 
tions which have discontinuities [5]. The latter is, unfortunately, the case in electromagnetic 
scattering problems. Suppose there is an interface between two media. It can be deduced 
from the Maxwell's equations that the components of the electric and magnetic fields, E t 
and Hj, tangential to the interface must be continuous, provided there is no surface electric 
current on the interface. From the Gauss law it follows that the components of the induc- 
tions, D t and B t , normal to the interface must be continuous, provided there is no surface 
charge on the interface. In contrast, the normal components of the fields and the tangential 
components of the inductions can be discontinuous. Their discontinuities are proportional to 
discontinuities of medium parameters (e.g., discontinuities in plasma frequencies in Lorentz 
models). Therefore, in either the induction or field representation, there are components 
which suffer discontinuities at the interface. Consequently, the convergence rate of the split 
algorithm for Maxwell's theory might be slower than that in quantum mechanics with a 
discontinuous potential because in the latter case the wave function remains continuous. 

Another source of errors that affects the convergence rate as N increases is the aliasing 
problem in the fast Fourier transform on the grid. Note that, even though the initial wave 
packet is band-width limited and the grid is chosen fine enough to eliminate errors in doing 
its fast Fourier transform back and forth, the wave packet looses this property after the 
operator exp(AtV Q ) is applied to it. As a result, the aliasing problem arises in spatial 
domains where V Q varies (typically at interfaces between different types of media). 

The above two problems that also reduce the accuracy of the algorithm are well known 
and studied in the theory of the fast Fourier transform [5]. The only way to cope with 
them is to make the grid finer in the areas where medium parameters have discontinuities. 
However, the fast Fourier algorithm requires a uniform equispaced lattice, which might lead 
to wasting computer resources if the increased resolution is necessary only in relatively small 
areas of the computational volume of the problem (e.g., only near an interface between two 
media). There are several ways to modify the algorithm when the above problems are too 
expensive to overcome by making a uniform grid finer. 

First, the grid can be made fine enough so that the action of powers of the Hamiltonian 
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TiP on the state vector ty® is sufficiently accurate in the Fourier basis as specified by the 
rules (3.4) and (3.5). The operator H Q is projected onto the Krylov space spanned by 
vectors (H Q ) k ^Q , k = 0, l,...,n. Its exponent (the evolution operator) is then computed 
by diagonalizing H Q instead of using the Lie- Trotter formula. Usually, it is sufficient to 
take a low dimensional Krylov space. This method is known as the Lanczos method [26]. A 
detailed study of the Krylov-Lanczos method as well as other similar pseudospectral methods 
in Maxwell theory will be given elsewhere. 

Second, one can give up a uniform grid, while preserving basic advantages of pseudospec- 
tral methods such as, e.g., exponential convergence. A possible way to emulate a non-uniform 
grid in a multiscale problem is to use wavelet bases. The problem here is to compute the 
action of exp(At7i^) in the split algorithm because the derivative operator V is not diagonal 
in this basis (in contrast to the Fourier basis). However, is expected to be sparse in a 
wavelet basis so that its direct diagonalization might not be expensive, and a significant 
reduction of computational costs can be achieved in the split algorithm, by using the fast 
wavelet transform, as compared to that in the Fourier basis. Otherwise, the use of (4.10) 
might be helpful in place of the direct diagonalization method. This approach has proved to 
be successful in solving multiscale initial value problems for the Schrodinger equation [27]. 
In the framework of Maxwell's theory for passive media, additional studies of several issues 
in time domain wavelet based algorithms, like, e.g., stability, would still be needed. 

Third, the fast Fourier transform algorithm remains in place but is applied to an auxiliary 
uniform grid that is related to a non-uniform grid in physical coordinates by a change of 
variables. Consider a change of variables y = y(x). A uniform grid in the new variables 
y would generate a non-uniform grid in the original Euclidean (physical) coordinates x. A 
desired local density of grid points in the physical space, to enhance the sampling efficiency 
in designated regions, can be achieved by an appropriate choice of the functions y(x) [28]. By 
necessity, the auxiliary grid spans a rectangular box (with periodic boundary conditions). Its 
pre-image in the physical space would not be a box in general, save for the case when the map 
y(x) splits into three individual one-dimensional maps yj = yj(xj), j = 1,2,3. Since, the 
derivatives are transformed as V x = -A(y) V y where the 3x3 matrix A is position dependent, 
the operator cannot be kept in the exponential. The action of its exponential on the 
state vector can be approximated by the leapfrog method in which only the action of on 
yjj^ is required. The latter can be done by the fast Fourier method according to the rules 
(3.5) and (3.4) applied to an operator being a product of position and derivative dependent 
operators. In contrast to the well studied quantum mechanical case, the algorithm appears 
to be unstable for media with absorption. In Section 9 a modification of the leapfrog scheme 
is proposed to achieve conditional stability. 
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9 The temporal leapfrog scheme 



Here we discuss a temporal finite difference scheme applied to the Maxwell theory for passive 
media in the Hamiltonian formalism. As has been pointed out, such a scheme might be 
helpful for reducing computational costs by using non-uniform grids in combination with 
some pseudospectral methods (e.g., wavelet bases or the fast Fourier method with a change 
of variables). A temporal finite difference scheme can be obtained by the following procedure. 
Let us integrate (1.1) over the interval (t,t + nAt). We have 

t+nAt / — 1 \ 

^t+uAt = ^ t + n J dr^ T =% + Atn (^2 c£ n) *t+kA t J + o(Ar +1 ) , (9.1) 

where the coefficients Cj^ used to approximate the integral are well known for any n as 
well as the accuracy of the approximation. For example, one can use the 3/8 Simpson rule 
for n = 3 or Bode's rule for n — 4. The iterating scheme allows one to compute the wave 
function at the sequential moment of time if it is known for n preceding moments of time. 
Only the simplest case n = 2, for which the mid-point approximation for the integral is 
taken, leading to Cf 1 = and cf^ = 2, will be considered in detail. It is also known as the 
leapfrog scheme: 

^t+At = %-At + 2AtH% • (9.2) 
The action of the Hamiltonian is computed in a suitable basis (as has been noted above). 
Apart from violation of the dispersion relation of electromagnetic waves, temporal finite 
difference schemes would generally be unstable in media with absorption, in contrast to the 
quantum mechanical case. The reason is that the Hamiltonian in (9.2) is not ant i- Hermit ian. 
Consequently, convergence to the continuum solution would also be violated. 

A general solution to (9.2) can be written in the form 

*nA t = (S^)" (9 ' 3) 

0$ = HAt ± Vl + H 2 At 2 , (9.4) 

for some initial state vectors and \l>At ( the vectors ^± are determined by them). Stability 
requires that there exists a positive constant C such that 

|*nAt| <C(|*+| + |*_|) , (9.5) 

for all < n < N, T = NAt and < At < r. Note that in general a solution of (1.1) 
may have a legitimate exponential growth if the hermitian part of the Hamiltonian, Ti + Ti,*, 
is not negative semidefinite (see Lemma 7.2). For this reason, a typical stability criterion 
would be equivalent to the condition [16] that there exists some positive constant K\ such 
that ll^^ll < 1 + K\At uniformly for all parameters of Q^ t and for < At < r, which is 
clearly the case for (9.4) if H is bounded. The latter leads to 



(G& ) ) n \\<e KlT 0.6) 
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so that a legitimate exponential growth of the solution is allowed, i.e., C ~ e KlT in (9.5). 
For passive media, the Hamiltonian satisfies the conditions of Lemma 7.2 and, hence, no 
legitimate exponential growth should be present in a numerical solution in order to achieve 
convergence. However, as we shall see shortly, the scheme (9.2) always generates an expo- 
nentially growing solution for media with attenuation. 

Let complex numbers z = Re t(f be eigenvalues of HAt. Since the spectral radius of 
does not exceed its norm, the necessary (von Neumann) condition to suppress an exponential 
growth of the solution (9.3) reads 



p{Q { ^) = max z ± Vl + z 2 < 1 . (9.7) 

The aim is to analyze the domain D of the complex plane for which (9.7) holds. Let 
f] = yl + z 2 and = r. The two inequalities in (9.7) require that for z G D, i? 2 +r 2 ±£ < 1, 
where £ = zrj+zf]. By combining the latter inequalities, one gets R 2 +r 2 < 1 or r 4 < (1— R 2 ) 2 . 
On the other hand, r 4 = 1 + R A + 2R 2 cos(2<^). Hence, cos(2</?) < —1 which is only possible 
if if = ±7r/2 or z = ±iR. The necessary condition (9.7) is satisfied if 

if = ±tt/2 , R 2 < 1 . (9.8) 

This does not yet guarantee that there is no norm growth. A norm growth, which is poly- 
nomial in time, can still occur. 

Let us investigate general properties of the solution of (1.1) when the Hamiltonian satisfies 
the von Neumann condition (9.8). For any matrix Ti there exists a similarity transformation 
so that S^TiS has the Jordan normal form. Let h z be a block of the Jordan normal form 
corresponding to an eigenvalue z of 7i. Any block hk is a q z x q z bi-diagonal matrix, q z > 1, 
with all the elements of the diagonal being equal to z and all the elements on the upper 
superdiagonal being equal to one. For q z — 1, h z — z is just a complex number. The norm 
of any solution of (1.1) cannot grow faster than || exp(t7i)||. Let a g^-dimensional vector <p t 
satisfy the equation d t (pt = h z <f> t . For a generic initial condition, the solution norm grows 
polynomially, \cf>t\ = 0(t 9z_1 ) as t — > oo. Using the similarity transformation S, we define 
the corresponding /z-norm of state vectors and the corresponding matrix norm (cf. (6.5)) by 
setting \i = S~ 1 *S~ 1 . From the equivalence of the norms || ■ || and || ■ \\^ (sec the proof of 
Theorem 7.3), the norm growth cannot be faster than 



\e tn \\ = \\e ts lns \\ = max||e t/lz || = 0(t 9_1 ) , q = m&xq z , t -> oo , (9.9) 

2 Z 



provided z = ±iR. However, a state vector norm growing polynomially with time is unac- 
ceptable from the physical point of view because in any passive medium there is no physical 
mechanism for such amplification of the field amplitudes in the large time limit. Conse- 
quently, we demand that any model Hamiltonian for a passive medium should be similar to 
a diagonal matrix (i.e., TC is diagonalizable) . In this latter case, the blocks h z of the Jordan 
normal form of H are just complex numbers z. Hence || exp(th z )\\ = \exp(±itR)\ = 1 so 
that the /x-norm of any solution of (1.1) is conserved according to (9.9). 
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Two important conclusions about the leapfrog scheme (9.2) follow from our analysis. 
First, the von Neumann condition (9.8) is also sufficient for stability. Indeed, if (9.8) holds 
then = WSg^S^W = p{G ( tt) = 1 and, hence, \\{G ( A t) n \\^ < 1 uniformly in n > 0. 

By the norm equivalence, ||(^^) n || is also bounded uniformly in n > 0. Second, reversing 
the argument, we conclude from the norm conservation of the stable leapfrog solution that 
no attenuation can be added to the Hamiltonian without destroying the stability and, con- 
sequently, the convergence to the continuum solution. Whenever the attenuation is added, 
the leapfrog solution would always contain an exponentially growing component, while this 
would not be so for a continuum solution (see Lemma 7.2). 

Since G^tG^t = 1> on ly one °f the two independent solutions in (9.3) would grow expo- 
nentially whenever the attenuation is added. Theoretically, for Ti + Ti* < the exponentially 
growing solution can be eliminated by choosing the initial condition so that = which 
is equivalent to the initial condition ty At = G^t^o- Practically, this is never possible due 
to rounding errors and/or numerical errors in computing G A + t^o- Even for a small \^-\ in 
(9.3) the growing part would eventually become comparable with the exponentially atten- 
uating solution generated by A reduction of the time step would not be helpful since 
the constant K\ in (9.6) is independent of At while the simulation time T is fixed by the 
dimension of the simulation volume and the speed of light. One needs at least to modify the 
scheme so that there exists a constant K p such that 

\\g^\\<l + K p At p , p>l, (9.10) 

for < At < r. Indeed, it follows from (9.10) that UG^TW < exp^TA**" 1 ) = 1 + 
0(At p_1 ) for all < n < N where iVAt = T. The norm growth could be reduced as much 
as desired by making the time step smaller. Next we show how to modify the leapfrog scheme 
to make (9.10) valid for at least p = 3 and, if the Hamiltonian is normal, an even stronger 
result holds, namely, K p = 0. 

Let H = Ho + V where V* + V < (negative semidefinite) and H* = -H . In (1.1) we 
make a substitution ^ t = exp(tV)$*. The new state vector $ t satisfies an equation with a 
time dependent Hamiltonian, 

dt$ t = e- tv H e tv $ t = Ht$t , (9.11) 

and with the same initial condition <3> = ^o- Applying the leapfrog method to (9.11) we get 
&t+At = ®t-At + 2At7-^$ t valid up to 0(At 3 ). Returning to the initial variables, we arrive 
at the following recurrence relation 

^t+At = C2AtVt-At + 2AtC At H ^t , (9. 12) 

where C At = exp(AtV). All the derivative operators are included into the anti-Hermitian 
part T~Lq of the Hamiltonian Ti, while the attenuation operator V might even be independent 
of position and, hence, C At has to be computed only once for given medium parameters and 



32 



time step. It can often be done analytically as, for example, in multiresonant Lorentz models 
(see Section 10). On the other hand, by Lemma 7.2, ||£At|| < 1 for any At > 0, and one 
might hope to stabilize the leapfrog scheme by satisfying the stability condition (9.8) for H 
only, that is, 1 + H^At 2 is positive semidefinite. This is indeed the case. The amplification 
matrix, ^t+At = GAt^t, for the recurrence (9.12), satisfies the equation 

Qai = C 2A tG A l + 2AtC At H . (9. 13) 

According to our analysis of the von Neumann stability condition (9.8), the anti-Hetmiticity 
condition of Ho in (9.12) and (9.13) can be weakened by demanding that Ho is related to 
an ant i- Hermit ian matrix by a similarity transformation. Some important properties of the 
amplification matrix obtained from (9.13) are stated in the following theorem. 

Theorem 9.1. Suppose there exists a similarity transformation such that S~ X HS = 
7~Ls+Vs where H* s = —Hs, the Hermitian part of Vs is negative semidefinite, V$+Vs < 0, and 
Hs also satisfies the von Neumann stability condition for the leapfrog scheme, 1 + H 2 s At 2 > 
(positive semidefinite). Consider the amplification matrix Q At of the modified leapfrog 
scheme (9.13) with Ho = SHsS^ 1 and V = SVsS^ 1 . Then there exists a norm || • || M 
equivalent to || • || such that Q At has the following properties: 

(A) if [W o ,V] = 0, 

WGlth < 1 (9-14) 

uniformly in n > 0; 

(B) if [Ho, V] 7^ 0, there exists a non-negative constant K 3 such that 

||0A*|l M < l + K 3 At 3 , (9.15) 

for < At < t and some positive r. 

Proof. Part (A). If Ho and V commute, the amplification matrix Q At = ^AtQ%t satisfies 
(9.13), provided Q\ t satisfies the same equation for V = (or C A t = 1), which one can 
easily check by substituting the solution into (9.13). Consider the norm associated with the 
similarity transformation S of the Hamiltonian, \\A\\fj, = WS^ASW- The norms || • || M and || • || 
are equivalent (see the proof of Theorem 7.3). Since Ho satisfies the von Neumann stability 
condition and is anti-Hermitian relative to the fi scalar product, = P(&At) = 1 

(according to the analysis after (9.7)). By Lemma 7.2, ||£AtlU < 1 ; and we infer that 
II^AtlU = Wi^tQltTh < H£a*||£ < 1 uniformly for n > 0. 

Part (B). Solving (9.13) by the perturbation theory in At, it is not hard to find that 

Qai - G At = At 3 JC At , Q\ t = C At / 2 G At C At / 2 , (9.16) 

where JC At is regular in the vicinity of At = and vanishes whenever Ho and V commute. On 
the grid, Ho and V are bounded operators. Hence we can find a constant K s = sup At ||/CaJ m 
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for some open interval < At < r. Making use of the inequality ||<?XJU — II^At^Hu < 1, 
we find 

II^AtlU = WGlt + At 3 /C At || M < 1 + At 3 K 3 , (9.17) 
which completes the proof. 

The norm deviation of the solution generated by the modified leapfrog scheme (9.12) 
from the stable solution generated by Q\ t is of order 0(At 2 ) for the entire simulation time 
T and, hence, by reducing At a possible norm growth can be suppressed as much as desired. 
Indeed, 

n-l 

\\nn i r>V \n\\ n \ rtn—kln nV \(nV \k\\ 

||t/At - (t/AtJ || M = \\ 2_^yAt yyAt-yAt){yAt) h 

k=i 

n—l 

< At 3 K \\Gm% < K 3 TAt 2 e K3 ™ 2 = 0(At 2 ) . 
k=i 

Since in the continuum limit At — > 0, both the amplification matrices Q\\ and Q^ t generate 
the same solution and all the powers of the former are uniformly bounded by construction, 
a natural question to ask is whether one can find a recurrence relation for the function 
^nAt = (^At) n ^o which could be used in place of (9.12). It is not difficult to derive an 
equation for Q\ t similar to (9.13), but, unfortunately, this equation cannot be converted into 
a simple recurrence relation for the wave function itself, like (9.12), suitable for numerical 
applications. 

It should be noted that if the operator C^t in the modified leapfrog scheme (9.12) is 
replaced by another £ s At such that £a< — £a< = 0(At 3 ) and — 1' then the convergence 

is not violated because Part B of Theorem 9.1 still holds. This observation is useful for 
analytic computation of C^ t . For example, in the conditions of Theorem 9.1, put S = 1. 
Let V = Vi + V2 so that both Vi,2 have their hermitian parts negative semidefinite. By using 
the split (4.2) we get 

C At = e AtV = e A t V 1 /2 e AtV 2e A t V 1 /2 + 0(At 3 } = £ ^ + (g^g) 

By Lemma 7.2, ||£ A J < 1 for At > 0. The operators Vi,2 can be chosen so that their 
exponentials can be computed analytically. 



10 Examples of the temporal leapfrog algorithm 

There are many possibilities to split the original Hamiltonian 7i into two parts that satisfy the 
conditions of Theorem 9.1 and thereby to make the leapfrog scheme stable and convergent. 
Basic guide lines for doing that are as follows. The Hamiltonian 7i should contain all 
the derivative operators in H and, yet, the von Neumann condition is easy to establish for 
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H, . It would also be helpful to have an analytic expression for C^t a t least up to order 
0(At 3 ). As an illustration, we discuss multiresonant Lorentz models and geometric optics. 
To distinguish between the splits of the Hamiltonian in the split and leapfrog algorithms, 
we shall use an index I ( "leapfrog" ) in the latter. 



I Lorentz models 

In the field representation of the Hamiltonian for multiresonant Lorentz models, we make 
the following decomposition 

Thanks to (2.28) and Ti^ = —Ho, the operator Ti^ is anti-Hermitian. From (2.29) it follows 
that the Hermitian part of Vf is negative semidefinite > 0). The exponential of Vf 
is easily computed according to (5.14) and (5.15). Let denote a six-component column 
whose three upper components coincide with ^ 2a ' 1 and three lower components equal £ 2a 
(see (2.9)-(2.11)). As a result we arrive at the following scheme 

<At = tf-At + 2Atn ^ F + 2 At ^V FMa C, (10.2) 

a 

% +At = e 2Atw - $_ M + 2At e At *- v MFa . (i . 3 ) 

Stability is ensured if Hf satisfies the von Neumann condition (9.8). Eigenvalues of Hf 
satisfy the equation 

det(z - H&) = z q det(z 2 - zH - V fm V M f) = (10.4) 

where the non-negative integer q depends on the number of resonances in the Lorentz model. 
Non-zero eigenvalues satisfy the so-called pencil equation whose theory is well developed 
and might be useful for more general models [29]. Here we shall find a simpler (practical) 
criterion sufficient for (9.8) to hold. Since the plasma frequencies may depend on position, 
we apply the following general idea [30]. Suppose we have a finite difference scheme with 
variable coefficients in space. Consider a corresponding finite difference scheme with frozen 
coefficients. It is obtained from the original scheme by fixing the coefficients to particular 
values everywhere in space. A finite difference scheme with variable coefficients is stable 
if all the corresponding finite difference schemes with frozen coefficients are stable [30, 16]. 
So let us fix the plasma frequencies to particular values. The spatial dependence of the 
eigenfunctions for the pencil problem in (10.4) is given by a harmonic factor exp(ik • x) and 
the corresponding eigenvalues are z = ±iy / c 3 k 2 "+~u^, where u) 2 p is defined in (5.16). Let 
Knax be the maximal norm of all wave vectors of the initial wave packet and Up iax be the 
maximal value of u p as a function of position, then a sufficient criterion for stability reads 

^c 2 k 2 max + < 1 . (10.5) 
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The scheme (9.12) becomes especially simple in the case of small attenuation, 7 a < uo a . In 
the complex representation of the auxiliary fields (2.31) (cf. (2.12)) the matter Hamiltonians 
H-Ma are diagonal and the action of its exponential is reduced to multiplication by a complex 
number e lUaAt (see Section 5). 

The stability condition (10.5) can be improved if one uses the induction representation 
arranging the split according to (4.11), that is, H, 1 ^ = Hq and V[ = V 1 . In this case the 
conditions of Theorem 9.1 are met if instead of (10.5) we demand a weaker condition 

Atck max < 1 . (10.6) 

To prove this, we note first that by the similarity transformation defined in (2.23) we get 

S-^(*° °) + ( v ° p V %)=K + * a . (10.7) 

Then (10.6) is obviously the von Neumann stability condition for 7i s , while the Hermitian 
part of V| is negative semidefinite if 7 a > 0. The scheme is obtained from (9.12) by replacing 
% -> *t, H -> and V V' as defined in (4.11). Since the left hand side of (10.7) 
coincides with H F , it can also be viewed as the leapfrog scheme in the field representation 
but with the split different from (10.1). This illustrates the point that the stability condition 
of the scheme (9.12) depends strongly on the choice of Tioi- The price for a simpler stability 
condition in the induction representation is the lack of an explicit form of £a<- However, 
this problem can be circumvented by making use of (9.18). Indeed, V| = V F + Vf where 
V F is defined in (4.7). The exponentials of these operators are computed in Section 5. We 
also have || exp(AtV F )|| = 1 because V F * = -V F and || exp(AtVf )|| < 1 by Lemma 7.2 for 
a non-negative At. We set 

so that H^AtH/i = II^^At^ll — 1- The operator (10.8) differs from C^ t = exp(AtV 7 ) 
= S cxp(AtVj.)5 _1 by terms of order 0(At 3 ) and, hence, according to (9.18), can be used 
in place of £a* in the leapfrog scheme without destroying its convergence and stability. 



II Geometric optics 

Another simple example is the case of geometric optics. For sake of simplicity we assume the 
medium to have no magnetic properties. A generalization is straightforward. Let e = e(x) 
be the dielectric constant of the medium. If the medium is not isotropic, then e is symmetric 
positive definite 3x3 matrix everywhere in space. We rewrite Maxwell's equations in the 
form 

= n G ^ t , H G = ( _ cV x ° (£ _ 1 } C ^ X ) , (10.9) 



36 



where the parentheses in (e *) mean that the induction is first multiplied by e 1 and then 
the curl of the resulting vector field is computed. Consider the scalar product 

tylti) = f dr , ^=( £ q 1 J)- ( 1(U0 ) 

In the grid representation of Section 3, the integral is replaced by the sum over grid points 
and Tic becomes a finite matrix. The Hamiltonian is anti-Hermitian with respect to this 
scalar product, Tt G /j, = —\iHg- Therefore the corresponding /j, norm is preserved in the time 
evolution generated by exp(tTi,G), that is, (V'/jV't) = (V'ojV'o)- The electromagnetic energy 
of the wave packet is conserved because it is proportional to the \i norm of the initial state 
vector. Consequently, we expect that for a sufficiently small At the original leapfrog scheme 
(9.2), 

l&At = ti-M + 2AfWG# , (10.11) 

becomes stable. To find a sufficient condition for stability, the same idea of finite difference 
schemes with frozen coefficients can be used. It obviously leads to a condition similar to 
(10.6), 

where k e max is the maximal norm of all wave vectors in the medium which can be estimated 
by a/ p{s)k max with k max being the maximal wave vector of the initial pulse in vacuum. The 
spectral radius p(e) is understood as the maximal spectral radius of e(x) over x. If the 
Fourier basis is used to compute the derivatives, the algorithm does not violate the Gauss 
law. However, the algorithm would not conserve the fi norm (or energy), rather a quantity 
which, in many cases, approximates the energy. Multiplying (10.11) by ty\ using the scalar 
product (10.10), we infer that 

«a*, ti) = ti-M) = ■■■ = O ■ (10.12) 

By expanding the exponential in ^ +At = exp(AiWc)^/ into a Taylor series in both sides of 
(10.12) and making use of the anti-Hermiticity of Hq, we find that the energy conservation 
violation is of order 0(At 2 ). Thus, it can be made as small as desired by reducing the time 
step. 



11 Conclusions 

The initial value problem in Maxwell theory for passive media has been reformulated in the 
Hamiltonian formalism. The path integral representation of the fundamental solution of the 
Hamiltonian evolution equation has been used to develop a time domain numerical algo- 
rithm for solving the initial value problem. The algorithm exhibits the main advantages of 
pseudospectral methods for solving differential equations such as an exponential convergence 
(and, hence, a greater accuracy), the absence of dispersive errors and numerical efficiency In 
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addition, the algorithm is unitary, meaning that the energy of the initial pulse is conserved 
whenever the medium attenuation vanishes (Theorem 6.1). For widely used multiresonant 
Lorentz models, the algorithm is unconditionally stable (Theorems 7.1 and 7.3), and, for a 
generic passive medium, conditional stability can always be achieved (Theorem 7.4). As the 
time step At goes to zero, the algorithm accuracy is of order 0(At 2 ) (Theorem 8.1). It is 
possible to increase the convergence rate (accuracy) up to any desired order 0(At n ), n > 2. 
However, computational costs for increasing the accuracy in such a way are not necessar- 
ily lower than those for decreasing the time step in the original algorithm. An important 
advantage of the algorithm is that the Gauss law holds exactly in the process of numerical 
simulations with no extra computational cost (Theorem 8.2). 

A drawback of the algorithm is related to well known problems of the fast Fourier method. 
Namely, a slower rate of convergence for non-smooth functions and aliasing. This might, 
perhaps, limit the advantages of the algorithm in some type of scattering problems with com- 
plex target geometries. Numerical tests are needed for a quantitative conclusion. There are 
several pseudospectral methods for approximating the fundamental solution of the Hamil- 
tonian evolution equation that can help to circumvent this problem. We have analyzed one 
of them and formulated its stability criteria in the case of general passive media (Theorem 
9.1). Numerical tests of the modified leapfrog scheme are presented in [4]. The results are 
compared with known theoretical and experimental studies of the system investigated (ex- 
traordinary transmission gratings [31]). Other methods will be discussed elsewhere as well 
as the case when radiation sources (antennas) are included. 

It is believed that the proposed algorithm would be useful in numerical studies of elec- 
tromagnetic pulse propagation in passive media (e.g., foliage, soil, etc), photonic crystals 
and devices, nonostructured materials, and also in applications to scattering problems with 
targets made of dispersive materials. 
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12 Appendix 



I Initial pulse configurations 

In principle, any field configuration can serve as the initial configuration. However it is 
often desired to have an initial pulse with some specific properties (bandwidth, polariza- 
tion, direction of propagation, etc). Yet, the initial wave packet should be built of radiation 
(propagating) electromagnetic fields. A simple method based on the fast Fourier trans- 
form algorithm is given below to obtain initial configurations made of radiation fields with 
designated properties. 

A general solution of the Maxwell equations in vacuum can be written in the form 

E t (r) = J dk (C +) (k) e * kr + icfe * + C ~ } (k) e lk ' r - lcfct ) = C[ +) (r) + C t (_) (r) (12.1) 

B t (r) = -jL= Vx (c< +) (r)-C<->(r)) . (12.2) 

The representation (12.2) for the magnetic induction follows from the Maxwell equations 
and that the complex amplitudes C^(k) satisfy the transversality and reality conditions 
which are, respectively, 

k-Cf)(k)=0, C^(k)=cW(-k). 

The representation (12.1) and (12.2) holds for any moment of time. Hence, we can set t — 
to generate suitable initial conditions for an electromagnetic pulse propagating in empty 
space by choosing specific functions C^^r). 

Consider a few examples. Let there be translational invariance along the y axis. In this 
case the fields depend only on x and z, i.e., r = (x,0,z). Accordingly, the wave vector has 
the form k = (k x , 0, k z ) and dk = dk x dk z . Let e 2 = (0, 1, 0) be the unit vector along the y 
axis. Introduce 

C±(r) = \e 2 A) e~^ 2 / 2 e^ k °' r = e 2 C*(r) , (12.3) 

where A and k are real constants, and k is a fixed wave vector. Note that the field Cq ± ' ) 
is automatically transversal. The corresponding fields determine suitable initial conditions 
to generate a pulse propagating in the direction of k whose frequency band is centered at 
ujq = ck and its width is proportional to ck. The pulse is linearly polarized along the y axis: 

E = e 2 E , E = C { +) + , (12.4) 

Bo = ^= V x e 2 (c { +) - C { ~ ] ) • (12.5) 

The action of the differential operator is defined via the fast Fourier transform (see Section 
3) on the grid fine enough to support the bandwidth limited function (12.3). In the Fourier 
basis, iV/^/^K -> -k/k. 
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To obtain suitable initial conditions for a pulse propagating in the direction k and whose 
polarization lies in the rrz-plane, we make use of the electromagnetic duality of Maxwell the- 
ory which states that the dynamics remains unchanged when electric and magnetic charges 
switch places and simultaneously E t — > — B t and B< — > E t . According to the duality theo- 
rem, we can take 



Finally, suitable initial conditions for a pulse propagating in the direction k with a 
generic elliptic polarization are obtained by taking a linear combination of the above two 
initial conditions for two independent linear polarizations of the pulse. The amplitudes 
Co^r) can also be set numerically from actual measurements of a particular pulse of interest. 

II Conductivity and absorbing boundary conditions 

In numerical simulations, the grid in coordinate space is of necessity finite. In scattering 
problems we are interested in the pulse shape and polarization which we wish to compute 
in the asymptotically large coordinate region. This requires that not only the leading edge 
of the reflected pulse should have reached the asymptotic region, but also the trailing edge 
should have done so as well. This is essential if the reflected pulse propagates in a highly 
dispersive medium, or the target has a complex shape, or both. A complication arises from 
the very nature of the fast Fourier transform method. The method is designed to describe 
periodic functions and, consequently, if the pulse has a finite amplitude at the edge of the 
grid, this finite value would appear back at the other edge, with totally disastrous results 
for the computation. In quantum computational physics this problem is often solved by 
using an optical potential that absorbs the signal as it reaches the grid boundary. A similar 
method can be developed for our treatment of Maxwell's theory. Before we do so let us 
point out that an absorbing boundary condition is not the only way to solve the problem. 
For instance, in the case of a complex target, an ancillary grid may be defined in one of the 
coordinates which extends to large distances. The pulse may be transferred in a gradual 
manner from the small grid (near the target) to this larger grid to prevent the pulse from 
ever reaching the edge of the small grid. This technique can also be applied to generate a 
pulse by an antenna of a complex construction. The dynamics of the portion of the pulse on 
the larger grid may be treated analytically (if dispersion properties of the medium are not 
too complex). 

In quantum mechanics absorbing boundary conditions are made by adding an imaginary 
potential to the Hamiltonian with support near the grid edges. In the Maxwell theory, the 
same can be achieved by adding conductivity which gradually increases as the grid edges are 




(12.7) 



(12.6) 
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approached. An interaction of conducting media with electromagnetic radiation is described 
by Ohm's law, 

J* = ciE t , (12.8) 

combined with Maxwell's equation (2.1), where the displacement current is amended as 
d t T>t — ► d t Dt + (47rcr/c)E t with a = er(r) being the conductivity of the medium. Consider a 
linearly polarized plane wave moving along the z axis. Let E u be the Fourier transform of the 
only component of the electric field E t . Disregarding for a moment any possible anomalous 
dispersion of the medium, we find that E w satisfies the equation 



d 2 z K{z) + 



o\z\ 



K{z) = . 



12.9) 



Equation (12.9) is identical to the stationary Schroedinger equation with an optical (absorb- 
ing) potential being proportional to cr(z). In simulations of quantum wave packets it has 
been found that one of the optimal potentials has the form [17] 



aiz 



(n + iy 



a„ 



{z/LY 



n > 2 , 



;i2.io) 



in the interval z G [0, L\ and a(z) = otherwise. So, our next task is to find an optimal 
constant a n such that that the conducting layer would not reflect or transmit electromagnetic 
energy in some designated frequency band. Maxwell's equations in a conducting medium 
are form-invariant under the scaling transformations 



OJ 



(12.11) 



where a and j3 are positive constants. If the conductivity <r(r) was found optimal for a 
frequency uj and over a length L, then the optimal conductivity for a frequency (5uj would 
be (3cr((3r) and the new length over which it is taken to act would be L/(3. 

Let us first study the reflectivity of the absorbing layer. Suppose, a(z) = cr 9(z) where 
9{z) is the Heaviside function. If a monochromatic linearly polarized wave coming from the 
negative z region has an amplitude one, then the reflected wave has the amplitude [19] 



1 



1 + 

The energy of the reflected wave is 



1 + 



47ri<7o 



1 +iq u 



R„= \^\ 2 = qll^ + 0{qt) 



(12.12) 



(12.13) 



R w increases as the ratio q w gets higher and is small if q£/4 « 1. Consider now a(z) which 
monotonically increases from z = in the positive z direction. Let = k u (z) be a local 
wave vector of the wave in the conducting medium, 



'12.14) 
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We shall argue that the reflection is negligibly small if the local wave vector does not signif- 
icantly changes over a distance of order k' 1 , that is, 

k w (z + 5z) - k w (z) 



« 1 , Sz ~ k- 1 . (12.15) 



Making a linear approximation in (12.15), we infer that 

\d z qUz)\ « (2c)- l u (l+ql(z)) 3/4 , (12.16) 

which must be valid for all values of q^z) including small ones when the reflection is 
small. Inequality (12.16) allows us to reverse the argument, that is, the reflection is small 
if \9 z <loj(z)\ « {2c)^ 1 uj. Let be an average conductivity over a layer of width L, 
Gl = L~ l J Q L dza(z). In particular, for (12.10), &l = a n . For a monotonically increas- 
ing function, the derivative can be approximated as |<9 2 g w (,2)| ps ctl/L. This leads to a 
necessary condition on conductivity to suppress the reflection, namely, 

Our analysis is valid if the higher derivatives of a(z) are not large. This condition requires 
that the exponent n in (12.10) should not be less than two to insure a smooth behavior at 

2 = 0. 

The transmission can be estimated as follows. Suppose the pulse occupies a compact 
region Q. Let £p be the pulse energy. The pulse looses its energy as it propagates through 
a conducting medium according to Ohm's law, so that 

c - l d t £p = -(2c)" 1 [ dr a E 2 t < -8na n £p , (12.18) 
Jn 

where an = maxj] a. Therefore the pulse energy decay can be bounded from above by 

Sp < e-^ tcn . (12.19) 

In the one dimensional case (12.10), an = ai/in + 1). For the time t = L/c needed 
for a pulse to get through the layer of width L, the attenuation should be large, that is, 
87tL<7l /c(n + 1) >£> 1. Thus, the necessary conditions to suppress both transmission and 
reflection (that is, to ensure an almost total absorption) of the pulse are 

(n + l)c u 2 L n 

<a L < - — . (12.20) 



8ttL 8ttc 

By changing the Hamiltonian H Q , the conducting layer can be included into the split 
or leapfrog algorithm. Since the conducting layer produces attenuation, the conductivity a 
must be included into the operator C&t m the modified leapfrog scheme. It is also possible 
to create an absorbing and non-reflecting layer by using a passive medium (e.g. a Lorentz 
model). The analysis of the medium properties would be similar to that for a conducting 
layer. In fact, using a layer of a passive medium would offer more flexibility in solving the 
grid boundary problem. 
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